NONPROFIT DATA

Create NPO Data Frame

# Import subset from NCCS BMF Data on nonprofits in Syr MSA

load( "../DATA/syr_orgs.Rda" )


# Subset - Extract poverty orgs

pov_orgs <- syr_orgs[which (  syr_orgs$NTEECC == "B61" | 
                              syr_orgs$NTEECC == "B82" |
                              syr_orgs$NTEECC == "B92" |
                            syr_orgs$nteeFinal1 == "E" | 
                            syr_orgs$nteeFinal1 == "F" |
                            syr_orgs$nteeFinal1 == "I" |
                              syr_orgs$NTEECC == "K30" |
                              syr_orgs$NTEECC == "K31" |
                              syr_orgs$NTEECC == "K34" |
                              syr_orgs$NTEECC == "K35" |
                              syr_orgs$NTEECC == "K36" |
                            syr_orgs$nteeFinal1 == "L" |  
                            syr_orgs$nteeFinal1 == "O" |
                            syr_orgs$nteeFinal1 == "P" |  
                            syr_orgs$nteeFinal1 == "R" |  
                            syr_orgs$nteeFinal1 == "S" |
                            syr_orgs$nteeFinal1 == "T"), ]

Geocode Poverty Organizations

print("Don't run me")

# download specific ggmap, register Google key

devtools::install_github("dkahle/ggmap")

register_google(key = 'AIzaSyAMrECdXqpOt443ms6nzh118uUsqC2lE6M',
                account_type = "premium", day_limit = 15000)

# subset orgs for gecoding address information

pov_orgs_gps <- subset(pov_orgs, select=c(EIN:INCOME))

pov_orgs_gps$whole_address <- do.call(paste, c(pov_orgs_gps[ 
                                      c("ADDRESS", "CITY","STATE",
                                        "ZIP5")], sep = ", ")) 

pov_orgs_gps <- subset(pov_orgs_gps, select=c(whole_address, EIN:INCOME))

for (i in 1:nrow(pov_orgs_gps)) {
  latlon = geocode(pov_orgs_gps[i,1]) 
  pov_orgs_gps$lon[i] = as.numeric(latlon[1])
  pov_orgs_gps$lat[i] = as.numeric(latlon[2])
}

# Make sure you find a way to then take any PO Boxes and place them in the 
# center of centroids, and only the ones that are NAs.

Drop missing values, create geojson file for poverty organizations

head(syr_orgs)
load("../DATA/pov_orgs_gps.Rda")
pov_orgs_gps_nona <- na.omit(pov_orgs_gps)

geojson_write( pov_orgs_gps_nona, geometry="point", file="../SHAPEFILES/pov_orgs.geojson" )
## <geojson>
##   Path:       ../SHAPEFILES/pov_orgs.geojson
##   From class: data.frame

Create Service Catchment Areas

# I am planning on cutting this section out in favor of doing the buffering a different way in the analysis section


class( pov_orgs_gps_nona )
## [1] "data.frame"
coordinates( pov_orgs_gps_nona ) = ~lat+lon

class( pov_orgs_gps_nona )
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
plot( pov_orgs_gps_nona )

# Create buffer zone around each point

pov_orgs_gps_nona_buff <- gBuffer( pov_orgs_gps_nona, width=.025, byid=TRUE )
plot( pov_orgs_gps_nona_buff, main="NPO service buffer zone (1.5 miles?)" )

class( pov_orgs_gps_nona_buff )
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

CENSUS DATA

Download & clean Census ACS 2011-2015 5-yr data

Census API Package

# Census API 
# install.packages("devtools")
# devtools::install_github("hrecht/censusapi")

censuskey <- "f8064a17832b439e690be95ab0e0ef9a78617584"

var.list <- c("NAME", 
              "B25077_001E", "B25003_001E", "B25003_002E", "B25003_003E",
              "B25002_001E", "B25002_002E", "B25002_003E", "B01003_001E",
              "B19013_001E", "B17020_002E", "B12006_001E", "B23025_005E",                      "B15003_017E", "B15003_018E",
              "B03001_001E", "B03002_003E", "B03002_004E", "B03002_005E",
              "B03002_006E", "B03002_007E", "B03002_008E", "B03002_009E", 
              "B03001_002E", "B03001_003E", 
              "GEOID" )

acs_2015 <- getCensus( name="acs5",
                       vintage = 2015,
                       key = censuskey,
                       vars=var.list, 
                       region="tract:*", 
                       regionin="state:36" 
                     )


# subset to Syracuse 
# sum(acs_2015$county == "067")
# sum(acs_2015$county == "075")
# sum(acs_2015$county == "053")



acs_2015_syr <- acs_2015[which (acs_2015$county == "067" | 
                                acs_2015$county == "075" |
                                acs_2015$county == "053"), ]

Rename Variables

# Race variables are single race, non-hispanic. For example, "white" is really "white non-hispanic." This is true for every race category listed except, obviously, hispanic.


rename( acs_2015_syr, 
       c("B25077_001E" = "mdn_hous_val", 
         "B25003_001E" = "tenure_tot", 
         "B25003_002E" = "tenure_own", 
         "B25003_003E" = "tenure_rent",
         "B25002_001E" = "tot_occup", 
         "B25002_002E" = "occupied",
         "B25002_003E" = "vacant", 
         "B01003_001E" = "tot_pop", 
         "B19013_001E" = "mhh_income",
         "B17020_002E" = "poverty",
         "B12006_001E" = "labor_part",
         "B23025_005E" = "unemploy",
         "B15003_017E" = "high_sch",
         "B15003_018E" = "ged",
         "B03001_001E" = "tot_race",
         "B03002_003E" = "white",             
         "B03002_004E" = "black", 
         "B03002_005E" = "am_ind",
         "B03002_006E" = "asian", 
         "B03002_007E" = "islander",
         "B03002_008E" = "other", 
         "B03002_009E" = "mixed",
         "B03001_002E" = "not_hispanic",
         "B03001_003E" = "hispanic"
      
           ))

acs_2015_syr <- rename.vars(acs_2015_syr,
            c("B25077_001E", "B25003_001E", "B25003_002E", "B25003_003E",
              "B25002_001E", "B25002_002E", "B25002_003E", "B01003_001E",
              "B19013_001E", "B17020_002E", "B12006_001E", "B23025_005E",                      "B15003_017E", "B15003_018E",
              "B03001_001E", "B03002_003E", "B03002_004E", "B03002_005E",
              "B03002_006E", "B03002_007E", "B03002_008E", "B03002_009E", 
              "B03001_002E", "B03001_003E"),  
            c("mdn_hous_val", "tenure_tot", "tenure_own", "tenure_rent",
              "tot_occup", "occupied", "vacant", "tot_pop", 
              "mhh_income", "poverty", "labor_part", "unemploy",
              "high_sch", "ged",
              "tot_race", "white", "black", "am_ind",
              "asian", "islander", "other", "mixed",
              "not_hispanic", "hispanic"))

# Create census race proportion variable

acs_2015_syr$porp_white <- (acs_2015_syr$white) / (acs_2015_syr$tot_race)

names( acs_2015_syr )

MAPS

Census Shapefiles

# Move this chunk to your shapefiles folder to show how you produced the geojson files

# I made another Rmd files with this chunk but also kept the original chunk here. Is that what you would like?

setwd( "../SHAPEFILES/" )

# Census tracts shapefiles of Syracuse MSA, changed to geojson files
Mad_county <- readOGR(dsn="Madison-ct", layer="tl_2010_36053_tract10")
## OGR data source with driver: ESRI Shapefile 
## Source: "Madison-ct", layer: "tl_2010_36053_tract10"
## with 16 features
## It has 12 fields
## Integer64 fields read as strings:  ALAND10 AWATER10
geojson_write( Mad_county, geometry="polygon", file="Mad_county.geojson" )
## <geojson>
##   Path:       Mad_county.geojson
##   From class: SpatialPolygonsDataFrame
class(Mad_county)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(Mad_county)
title(main = "Census Tracts in Madison County, NY")

names(Mad_county)
##  [1] "STATEFP10"  "COUNTYFP10" "TRACTCE10"  "GEOID10"    "NAME10"    
##  [6] "NAMELSAD10" "MTFCC10"    "FUNCSTAT10" "ALAND10"    "AWATER10"  
## [11] "INTPTLAT10" "INTPTLON10"
On_county <- readOGR(dsn="Onondaga-ct", layer="tl_2010_36067_tract10")
## OGR data source with driver: ESRI Shapefile 
## Source: "Onondaga-ct", layer: "tl_2010_36067_tract10"
## with 140 features
## It has 12 fields
## Integer64 fields read as strings:  ALAND10 AWATER10
geojson_write( On_county, geometry="polygon", file="On_county.geojson" )
## <geojson>
##   Path:       On_county.geojson
##   From class: SpatialPolygonsDataFrame
class(On_county)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(On_county)  
title(main = "Census Tracts in Onondaga County, NY")

names(On_county)
##  [1] "STATEFP10"  "COUNTYFP10" "TRACTCE10"  "GEOID10"    "NAME10"    
##  [6] "NAMELSAD10" "MTFCC10"    "FUNCSTAT10" "ALAND10"    "AWATER10"  
## [11] "INTPTLAT10" "INTPTLON10"
Osw_county <- readOGR(dsn="Oswego-ct", layer="tl_2010_36075_tract10")
## OGR data source with driver: ESRI Shapefile 
## Source: "Oswego-ct", layer: "tl_2010_36075_tract10"
## with 30 features
## It has 12 fields
## Integer64 fields read as strings:  ALAND10 AWATER10
geojson_write( Osw_county, geometry="polygon", file="Osw_county.geojson" )
## <geojson>
##   Path:       Osw_county.geojson
##   From class: SpatialPolygonsDataFrame
class(Osw_county)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(Osw_county)  
title(main = "Census Tracts in Oswego County, NY")

names(Osw_county)
##  [1] "STATEFP10"  "COUNTYFP10" "TRACTCE10"  "GEOID10"    "NAME10"    
##  [6] "NAMELSAD10" "MTFCC10"    "FUNCSTAT10" "ALAND10"    "AWATER10"  
## [11] "INTPTLAT10" "INTPTLON10"

Add Data to Shapefiles

Add Syracuse MSA census tract data to census tract shp files

# drop NAs
Mad_county <- na.omit(Mad_county)
On_county <- na.omit(On_county)
Osw_county <- na.omit(Osw_county)

#Merge spatial files 
syr_msa <- union(Mad_county, On_county)
syr_msa <- union(syr_msa, Osw_county)
plot( syr_msa )

# Merge data with MSA shp file
syr_merged <- geo_join(syr_msa, acs_2015_syr, "GEOID10", "GEOID")
plot( syr_merged )
title(main = "Syracuse MSA, NY")

# Create geojson file of Syr MSA 
geojson_write( syr_merged, geometry="polygon", file="../SHAPEFILES/syr_merged.geojson" )
## <geojson>
##   Path:       ../SHAPEFILES/syr_merged.geojson
##   From class: SpatialPolygonsDataFrame

Create centroid points for each census tract

syr_merged_cen = gCentroid(syr_merged,byid=TRUE)

geojson_write( syr_merged_cen, geometry="point", file="../SHAPEFILES/syr_merged_cen.geojson" )
## <geojson>
##   Path:       ../SHAPEFILES/syr_merged_cen.geojson
##   From class: SpatialPoints
plot( syr_merged )
points(coordinates(syr_merged),pch=20)
points(syr_merged_cen,pch=20)

Create zoom enabled map of Syracuse, NY

map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY")

map %>% addProviderTiles(providers$CartoDB)

Set the view of the Syracuse MSA

# syr_merged_cen <- geojsonio::geojson_read("syr_merged_cen.geojson", what = "sp")

syr_map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr_merged,
              fillOpacity = 0.3
              ) %>%
  addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat) %>%
  addCircles(data = syr_merged_cen)

syr_map %>% addProviderTiles(providers$CartoDB)

Edit color format of polygons and data points

syr_map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr_merged,
              fillOpacity = 0.3, weight = 2
              ) %>%
  addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat,
             color= "red", opacity = 1.0) %>%
  addCircles(data = syr_merged_cen, color= "black", opacity = 1.0)

syr_map %>% addProviderTiles(providers$CartoDB)

Adding census tract information to map

# By race

pal <- colorNumeric(
  palette = "Blues",
  domain = acs_2015_syr$porp_white )


syr_map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr_merged, fillOpacity = 0.3, weight = 2,
              color = ~pal(acs_2015_syr$porp_white)) %>%
  addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat,
             color= "red", opacity = 1.0) %>%
  addCircles(data = syr_merged_cen, color= "black", opacity = 1.0) 

syr_map %>% addProviderTiles(providers$CartoDB)
# By Income level

pal2 <- colorNumeric(
  palette = "Blues",
  domain = acs_2015_syr$mdn_hous_val )


syr_map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr_merged, fillOpacity = 0.3, weight = 2,
              fillColor = ~pal2(acs_2015_syr$mdn_hous_val)) %>%
  addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat,
             color= "red", opacity = 1.0) %>%
  addCircles(data = syr_merged_cen, color= "black", opacity = 1.0) 

syr_map %>% addProviderTiles(providers$CartoDB)

ANALYSIS

Spatial analysis using leaflet

# URL's below explain how to use functions and arguments in leaflet. It also contains info which notes that the radius of a circle is in meters. Therefore, 1.5 miles converts to 2400 meters.
# https://rstudio.github.io/leaflet/shapes.html
# https://github.com/Robinlovelace/Creating-maps-in-R/blob/master/intro-spatial.Rmd

# Creating a buffer zone, 1.5 miles from every NPO
# I followed the info from this site and tried to apply it to the analysis
# http://walkerke.github.io/2016/12/spatial-pipelines/


# Example using one point: 
test_buffer <- c(-76.16674, 43.03231) %>%
  st_point() %>%
  st_sfc(crs = 4326) %>%
  st_transform(32618) %>%
  st_buffer(dist = 2400)

test_buffer_map <- test_buffer %>%
  st_transform(4326) %>%
  as("Spatial")

syr_map %>%
  addPolygons(data = test_buffer_map) %>% 
  addProviderTiles(providers$CartoDB)
# Attempt to use every pov_org location, but doesn't work. Coding error.
# pov_orgs_buff <- c(pov_orgs_gps_nona$lon, pov_orgs_gps_nona$lat) %>%
  # st_point() %>%
  # st_sfc(crs = 4326) %>%
  # st_transform(32618) %>%
  # st_buffer(dist = 2400)

# pov_orgs_buff_map <- pov_orgs_buff %>%
  # st_transform(4326) %>%
  # as("Spatial")

# syr_map %>%
  # addPolygons(data = pov_orgs_buff_map)

# Alternative way to create buffer for pov_orgs

pov_orgs_buff <- gBuffer(spgeom = pov_orgs_gps_nona, width = 2400)

# proj4string(pov_orgs_buff) <- NA_character_
# proj4string(syr_merged_cen) <- NA_character_
# proj4string(pov_orgs_buff) <- CRS("+proj=longlat +datum=WGS84")
# proj4string(syr_merged_cen) <- CRS("+proj=longlat +datum=WGS84")

# inter <- gIntersection(pov_orgs_gps_nona, syr_merged_cen)

# Chloropleth Map of median housing values

pal2 <- colorNumeric(
  palette = "Blues",
  domain = acs_2015_syr$mdn_hous_val )

syr_map <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
  addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
  addPolygons(data = syr_merged, fillOpacity = 0.3, weight = 2,
              fillColor = ~pal2(acs_2015_syr$mdn_hous_val)) %>%
  addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat,
             color= "red", opacity = .3, radius = 2400 ) %>%
  addCircles(data = syr_merged_cen, color= "black", opacity = 1.0) 

syr_map %>% addProviderTiles(providers$CartoDB)
# Intersection with Census data, capture data
  # number of orgs within 1.5 miles of tract's center 
    #  nearest neighbor?
    # distance analysis 


# Poportion of census tract that are high/low diversity
  # number of nonprofits by type in high/low racially diverse areas
  # number of nonprofits by type in low/high income census tracts


# Corr analysis between nonprofit access and census tract attributes
  # Number of nonprofits location by census tract income or race
    # Race
    # Income

Final section: Account for all urban census tracts in 25 cities